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ABSTRACT 

The aim of the present paper is to investigate the spatial structure of a protoplanetary 
disc whose dynamics is governed by magnetorotational turbulence. We perform a 
series of local 3D chemo-radiative MHD simulations located at different radii of a 
disc w hich is twice as massive as the standard minimum mass solar nebula of |HayasTii| 



( 1981 ). The ionisation state of the disc is calculated by including coUisional ionisation, 
stellar X-rays, cosmic rays and the decay of radionuclides as ionisation sources, and 
by solving a simplified chemical network which includes the effect of the absorption 
of free charges by |j,m-sized dust grains. In the region where the ionisation is too low 
to assure good coupling between matter and magnetic fields, a non-turbulent central 
"dead zone" forms, which ranges approximately from a distance of 2 AU to 4 AU from 
the central star. The approach taken in the present work allows for the first time to 
derive the global spatial structure of a protoplanetary disc from a set of physically 
realistic numerical simulations. 

Key words: accretion, accretion discs - planetary systems:protoplanetary discs - 
turbulence - instabilities - magnetic fields - MHD - radiative transfer 



1 INTRODUCTION 

According to the present understanding, the process of 
planet formation takes place in protoplanetary discs, which 
are accretion discs around young, solar-type stars of the 
T Tauri class. Obtaining knowledge about the physical con- 
ditions in protoplanetary discs is therefore of essential im- 
portance for developing viable theories of how planets may 
form. 

One of the most intriguing fact about protoplanetary 
discs is their rather short lifetime of only about 10 million 
years ( Hartmann et al.|1998 l. The big question is: How can 
the matter in the disc get rid of its angular momentum in 
such a short time-scale? The most likely explanation for this 
is hydromagnetic turbulence initiated by the magnetorota- 
tional instability (MRI), a process which was introduced into 
the context of accretion disc physics by |Balbus fc Hawley] 
( 1991 1. Numerical simulations show that magnetorotational 



turbulence leads to fast outward angular momentum trans- 
port that can explain the high accretion rates observed for 
protoplanetary discs ( King et al.||2007 |. However, the MRI 
will work properly only if there is good coupling between 
matter and magnetic fields. Since protoplanetary discs are 
cool objects, with temperatures in the range of several hun- 
dred down to a few tens of Kelvin in most parts of the disc, 
this is a critical issue. 



Where in the disc the MRI is active, and where not, does 
depend on the value of the ionisation level there. While in 
the hot, inner regions (at a distance R < 1 AU from the 
central star) coUisional ionisation suffices to provide good 
coupling, this is no longer true in the cooler regions further 
away from the star. There, the disc has to rely on other ion- 
isation sources, like the decay of radionuclides, cosmic rays 
and X-rays, which are emitted from the corona of the cen- 
tral star. In the planet-forming region, at distances of several 
AU from the star, the density is so high that neither the X- 
rays nor the cosmic rays are able to reach the midplane of 
the disc, leading to a poorly ionised, non-turbulent central 
"dead zone" (Gammie 1996 Armitage 20111. At these in- 



termediate distances, only the upper layers of the disc are 
expected to be turbulent and to still provide a small amount 
of angular momentum transport. In the outer regions, the 
surface density is low enough for the X-rays to penetrate 
the whole disc column, so the midplane becomes turbulent 
again. 



While analytical calculations like in Gammie ( 1996 1 or 



simplified 1-l-lD models like the one developed in |Kretke fc| 
Lin ( 2010 1 can already provide useful models for the struc- 
ture of protoplanetary discs, the definite answers on the 
questions of the size of the dead zone and the strength of 
the angular momentum transport in the disc can only come 
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from numerical simulations. At the present stage it seems 
very difficult to perform global simulations including real- 
istic physics due to the large computational cost and the 
numerical complexity (see Fromang & Nelson||2006| |2009 



The physical equations are then given by 



[Dzyurkevich et al.||20l0l |F ock et al.||2011| for examples of 
global protoplanetary disc simulations). On the other hand, 
the simpler local simulations, which model only a small part 
of the disc, do already at the present time allow for the inclu- 
sion of additional physics like radiation transport and disc 
chemistry ( |Flaig et al.|2010[ [ffirose fc Turner|2011 l. 

In the present paper, we investigate the spatial struc- 
ture of a protoplanetary disc by performing a series of local 
3D magnetohydrodynamical simulations located at differ- 
ent radii, including both radiation transport and the effect 
of a finite Ohmic resistivity. We choose optimistic physical 
parameters in order to obtain a small dead zone that fits in- 
side the domain that is simulated. Using this method, we are 
able to obtain a comprehensive picture of both the vertical 
and the radial structure of a magnetorotationally turbulent 
protoplanetary disc. 

It should be noted that apart from Ohmic resistivity, 
ambipolar diffusion might also reduce the saturation level of 
the MRI. The strength of this effect depends on the value 



of the neutral-ion collision frequency (see, for example, Bai 
|fc Stone||2011[ ). In the fully turbulent regions of our model, 
the ratio of collision frequency to orbital frequency is > 100, 
suggesting that in these regions, the saturation level would 
not be strongly affected by ambipolar diffusion. 

The plan of our paper is as follows: In Sec. 2 we de- 
scribe our physical model and the numerical setup. Sec. 3 
presents the results of the numerical simulations and draws 
connections with astrophysical observations. In Sec. 4, we 
conclude. 
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with the total energy etot ~ p/ij— 1) + pv^ /2 + B'^ /2fio, the 
total pressure ptot = p + /2fio, and 

/cxt = — 2pJ7z X V + ipQ^x X — pfP z z (1) 

denotes the source terms arising in the local shearing box 
frame ( jFlaig et al.|2010"| , with fi the local orbital frequency. 
The radiation flux is given by F = -(Ac/Kp)Var*. The use 
of the flux-limiter A (for which we use the form suggested by 
[Levermore fc Pomraning||1981[ ) makes the radiation trans- 
port method applicable also to optically thin regions. 

The above equations are solved using a conservative fi- 
nite volume scheme. The scheme employs the HLLD Rie- 
mann solver of |Miyoshi fc Kusano| ( |2006[ ), which yields a 
high effective resolution at moderate computational cost. 

The value of the resistivity rj is calculated by including 
various ionisation source and by solving a simplified chemi- 



cal network. As in Hirose & Turner ( 2011 1, the values for the 



resistivity are read from a precomputed table. The value of 
the resistivity depends on the density, the temperature and 
the local ionisation rate due to X-rays, cosmic rays and the 
decay of radiormclides. We now describe the prescription ac- 
cording to which the resistivity is calculated. 



2 MODEL SETUP 

Our basic setup is very similar to that of the radiative pro- 



toplanetary disc simulations described in Flaig et al. (20101 



The simulations take place in the so-called stratified local 
shearing box, which is a rectangular box that covers the full 
vertical height of the disc but has only a small radial and 
azimuthal extent. This allows the use of local Cartesian co- 
ordinates {x,y,z), where x corresponds to the radial, y to 
the azimuthal and z to the vertical direction, respectively 
(for more information on the shearing box setup, sec Haw- 
ley et al.|p95} [Stone et al.1[T996l [Stone fc Ga~diner,,2010p . 
At the vertical boundaries, outflow boundary conditions are 
applied, that allow matter and radiation to escape from the 
disc (see [Flaig et al.|2010[ ). 

The disc gas is described by the equations of magneto- 
hydrodynamics, where we include radiation transport in the 
one-temperature flux-limited diffusion approximation ( [Flaig[ 
et al.|20l0 1 as well as the effect of a finite Ohmic resistivity. 



2.1 Ionisation sources 

2.1.1 Collisional ionisation 

In the hot, inner region inside 1 AU, collisional ionisation 
is the dominant ionisation source. The ionisation level aris- 
ing from this process is calculated using the Saha equation, 
which is given by 
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exp(-25 188/r) 
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(2) 



where rio and Un are the electron and neutral number den- 
sities in cm~^, respectively, T is the temperature given in 
Kelvin and we have assumed a potassium abundance of 10~^ 
(see, for example |Fromang et al.||2002| ). 



2.1.2 Stellar X-rays 

In the region where collisional ionisation is low, stellar X- 
rays are the dominant ionisation source. For the ionisation 
rate due to X-rays, we use the formula as given by [Turner| 
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where ^^(i?,^) are the column densities at radius R above 
and below a given point z, 

r'ltoo 



p{R, z') dz' , 



(4) 



and I/XR is the stellar X-ray luminosity which we take as 



2 ■ lO-" 



(5) 



Although this is a rather optimistic value, it is still inside the 
usual range of 10^' — 10"^^ ergs^^ for the X-ray luminosities 
observed for young stellar objects. 




Figure 1. Schematic diagram of tlie charge-transfer reactions in 
the ODD network. The symbols m, M and gr denote the gen- 
eralised molecule species (i.e. mainly H2), the generalised metal 
species (mainly K) and the dust grains, respectively. The red ar- 
rows denote the ionisation of molecules, the black arrows denote 
the recombination and charge transfer reactions involving only 
molecules and metals, while the blue arrows denote the absorb- 
tion of charged particles on the surface of dust grains. 



2.1.3 Cosmic rays 

Another possibly important ionisation source are cosmic 
rays, i.e. highly energetic particles from the interstellar 
medium that hit the disc. If these particles can reach the in- 



terior of the disc, their contribution is given by ( Umebayashi 
fc Nakano|[2009t 
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where ^cr,o is set to the cosmic ray ionisation rate in 
the interstellar medium, Ccr,o ~ 10~^^s~^, and xcr ~ 
96.0gcm-^ 



2.1.4 Decay of radionuclides 

Finally, we also include the ionisation arising from the decay 
of radionuclides. Long-lived radionuclides provide a back- 
ground ionisation level of about 10~^^s~^, where the dom- 
inant contribution is due to '^''K. This rate can be signifi- 
cantly increased if one includes the effect of short-lived ra- 
dionuclides, which are no more present in the solar system. 
The most important contribution comes from the decay of 
^®A1, which yields an ionisation rate of ("ra ~ 9.2- 10"^'' s~^ 
based on the mean interstellar abundance and eight times 
this value for the projected abundance of the young so- 
lar system, which is based on the ratio of different alu- 
minum isotopes found in CAIs in meteorites (Umebayashi, 
fc Nakano|2009 1. Here we adopt an ionisation rate of ("ra = 
7-10-^«s-\ 



2.2 Chemical network 



In the region where coUisional ionisation is ineffective, the 
ionisation state of the disc gas is determined by a balance 



between ionisation due to the ionisation sources discussed 
above and recombination of free charges inside the gas. We 
use the extended Oppenheimer-Dalgarno network of |Ilgner] 
I&: Nelso"n| ( |2006 !) , which approximates the gas-phase chem- 
istry by a generalized molecule species m and metal species 
M, extended by the introduction of spherical dust grains 



which have a grain mass density of Pg 



3 gcm"'^'. This 



extended network is called the "ODD network" in the fol- 
lowing. 

The charges of the newly introduced dust grains are 
tracked in the range from —2 to 2, leading to five new 
species gr'^" to gr'^"'". Furthermore, two mantle-species are 
introduced m(gr) and M(gr), which represent molecules and 
metals frozen out on the grain surfaces, leading to a total 
number of 12 species. 

The reactions among the different species can be found 
in ( Ilgner fc Nelson||2006 table 1, 3 and 4) and are graph- 
ically illustrated in Fig. [IJ (As we are mainly interested in 
the overall ionisation fraction, only the reactions involving 
charged particles are shown). The red ionisation reaction 
and the black charge transfer and recombination reactions 



are the same as in the Oppenheimer & Dalgarno ( 1974 1 net- 



work, whereas the reactions involving dust grains are drawn 
in blue. 

The resistivity is calculated from the ionisation fraction 
according to the usual formula 



77 = 234a;o/\/T. 



(7) 



The main effect of the dust grains is to provide new, indirect 
recombination paths for electrons: The grains sweep up free 
electrons, charging up negatively and simultaneously acquire 
positive charge through charge transfer reactions with metal 
and molecule ions. Grain-grain charge transfer reactions en- 
sure that oppositely charged grains neutralise themselves 
resulting in a relatively narrow grain charge distribution. In 
contrast to the metal ions, which provide a charge reser- 
voir that can only have a significant effect if enough metal 
is present in the gas, dust grains behave as a kind of cata- 
lyst for recombination. It is therefore possible that even tiny 
amounts of dust grains can change the equilibrium electron 
fraction significantly. 
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Figure 2. Equilibration time of the ODD network over orbital 
time as a function of the disc position. The contour lines corre- 
spond to a magnetic Reynolds number of 10^ and 10*, respec- 
tively. For further comments, see the text. 



2.2.1 Equilibration timescale 

Fig. [2] shows a plot of the equilibration timescale tequi over 
the orbital time tdy-n ~ 2n/0 for the ODD network, with a 
grain size of 10 jxm, a dust-to-gas ratio of lO"*, a metallicity 
of IQ-i" and a disc temperature that varies with the distance 
R to the central star according to 

T{R) = 280 (i?/AU)"'/^ K. (8) 

In addition, we also show contour lines for the magnetic 
Reynolds number defined as 

Re =M^-£l. 

77 r] U 

where Cg is the sound speed and H = 
scale height. The critical value of Ren 



(9) 



- Cs/O the pressure 
for which the tran- 
sition from the MRI-active to dead occurs, lies somewhere 
between 10^ (for the case where the vertical magnetic fiux 
through the disc is non-zero) and lO** (for a zero net flux con- 
figuration) (see Ilgner fc Nelson|2006 1 . As can be seen from 
Fig.[2j the equilibration timescale mostly stays below the or- 
bital time in the region where the transition from MRI-dead 
zone to the active zone occurs. This means that turbulent 
mixing effects can be neglected and the disc gas can be as- 
sumed to be in chemical equilibrium all the time. We note 
that when using a more complex chemical network, the equi- 
librium timescale is expected to be even shorter, since there 
are then more recombination channels available ('Bai"2011'). 
The chemical equilibrium value for the electron concentra- 
tion depends only on three parameters, the ionisation rate (, 
the temperature T and the total number density of the gas 
n (assuming a constant dust-to-gas mass ratio and metallic- 
ity). This fact allows us to use a previously computed, three 
dimensional look-up table for the equilibrium resistivity val- 
ues in our MHD simulation instead of evolving the chemical 
network during the simulation. 

2.3 Physical parameters 

We take the central star to have one solar mass. Concerning 
the surface density profile, we assume a disc mass which 



Parameter 


Symbol 


Value 


Mass of central star 


M* 


IM0 


Adiabatic index 


7 


1.4 


Mean molecular weight 




2.35 


X-Ray luminosity 




lO'^i ergs~l 


Cosmic ray ionisation rate 


CcR 


10-17 s-l 


Ionisation due to radionuclides 


Cra 


7- 10-i9s-i 


Dust grain size 


agr 


lOn-m 


Dust-to-gas ratio 


Pd/P 


10-* 


Material density of dust grains 


Pgr 


3 gcm"^ 


Metallicity 


M 


10-10 



Table 1. Summary of basic physical parameters of the physical 
model. For further comments see the text. 



is two times the disc mass as given by the minimum mass 
solar nebula model oflHayashi ( 1981 1, i.e. the surface density 



profile is given by the formula 



Eo = 3400 g cm^ 



R 
lAU 



-3/2 



(10) 



For this choice of surface density, the disc is to a large part 
optically thick in the whole radial range spanned by our sim- 
ulations (see Sec. 3). This means that the one-temperature 
approximation is well satisfied. 

The disc gas is characterised by an adiabatic index of 
1.4 and a mean molecular weight of 2.35. In accordance with 
astrophysical observations, we assume that the dust particles 
in the disc have grown to a size of several microns, choosing 
a monodisperse dust grain population with grain size Ogr — 
10 |J.m for the chemical network. We assume that the small 
dust grains are partly depleted due to grain growth and 
assume a dust-to-gas ratio of pd/p ~ 10-*. The metallicity 
is taken to be M = lO-i". 



As in [Flaig et aL] ( |2010[ ), we use the [Bell fc Lin| ( |1994 l 
opacity model. Since it is based on a dust grain population of 
smaller size than that used for the chemical network, these 
two choices are not fully consistent. This fact is, however, 
not very serious, since the opacity does depend only weakly 
on the grain size ( Pollack et al.|[T985 l and the opacities in 
protoplanetary discs are not well known anyway. For this 
reasons, the use of an improved opacity model is deferred to 
a future work. 

For reference, all the basic physical parameters used in 
our model are summarised in Table [TJ 



3 SIMULATION RESULTS 

We perform seven simulations located at different radii, 
starting from a distance of 1 AU from the central star up 
to a distance of 7AU. The simulations are initialised us- 
ing the same set of initial conditions as described in |Flaig| 
|et al.| (|2010[ ), with an initial constant temperature To (see 
Table [2|), a hydrostatic (Gaussian) density distribution with 
a surface density according to Eq. ( 10 1 and a zero net flux 



magnetic field that has an initial plasma beta of /3 = 100 at 
the midplane. In terms of the initial pressure scale-height, 
the box size is one scale-height in the radial direction and 
6 scale-heights in the azimuthal direction. The vertical box 
size is 14 scale-heights for the simulations at i? = 1 AU and 
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R [AU] 


Box size [AU] 

Lx X Ly X Lz 


Resolution 


Runtime [y] 


To [K] 


*oool 
torb 


H [AU] 


m) [K] 




(("» 


1 


0.07x0.41x0.97 


32x64x384 


100 


1200 


42.3 


0.07 


1537 


1.02 ■ 10^ 


0.038 


2 


0.11x0.67x1.35 


32x64x256 


283 


400 


29.8 


0.09 


249.7 


441.5 


9.3 ■ 10-5 


3 


0.15x0.88x1.17 


32x64x256 


519 


200 


24.0 


0.11 


139.0 


515.9 


0.001 


4 


0.22x1.35x1.80 


32x64x256 


800 


200 


6.5 


0.15 


122.5 


1091.8 


0.004 


5 


0.31x1.89x3.15 


32x64x256 


1118 


200 


2.3 


0.28 


161.8 


24903.4 


0.013 


6 


0.36x2.15x3.58 


32x64x256 


1469 


150 


6.7 


0.34 


163.1 


17336.4 


0.039 


7 


0.45x2.71x3.61 


32x64x256 


1851 


150 


3.3 


0.24 


54.4 


1983.14 


0.011 



Table 2. Overview of the simulation runs that were performed. Listed are the distance to the central star (first column), the box size 
(second column), the number of grid cells in each direction (third column), the time in years each simulations has run (fourth column), 
the initial temperature (fifth column), the initial cooling time in units of the orbital period (sixth column), the final scale height (seventh 
column), the mean temperature (eighth column), the mean magnetic Reynolds number (ninth column) and the alpha parameter (tenth 
column). The time-averages for the quantities listed in columns eight to ten have been performed from 40 to 100 orbits. 



10 scale-heights for the simulations at 7? = 5 and R = 6 
AU, while for all other simulations it is 8 scale-heights. The 
corresponding box sizes as measured in AU can be found 
in Table (2] The simulations are seeded with small random 
velocity perturbations of order ~ O.OlCs, leading to rapid 
development of the MRI, such that the disc reaches a fully 
turbulent state already within the first ten orbits (see Fig. [4] 
later in the paper). Each simulation ran for 100 local or- 
bits. The corresponding runtimes in years are also listed int 
Tabled 



3.1 Time history 

3.1.1 Thermal equilibrium 

In our setup, where we neglect the heating due to the irra- 
diation from the central star, the thermal structure of the 
disc is determined by a dynamical balance between internal 
heating and radiative cooling. Estimating the cooling time 
due to radiative diffusion from the relation 



1700 



tc 



(11) 



where the radiative diffusion coefficient Drad is given by 

i jFlaig et al.|20l6 l, we expect the simulations to reach ther- 
modynamical equilibrium after a few tens of orbits (see Ta- 
ble [2]| . The resulting temperature should then be indepen- 
dent of the initial temperature. 

In Fig. [3] the temporal evolution of the density- weighted 
spatially averaged temperature 



J pTdV 
JpdV 



(13) 



is shown. This plot suggest that the simulations do indeed 
reach a state of at least approximate thermodynamical equi- 
librium, with the possible exception of the simulations at 



^ Note that the heating due to the turbulent dissipation of ki- 
netic and magnetic energy is included in this, since by virtue of 
the conservative nature of the underlying numerical scheme, any 
loss of kinetic and magnetic energy is automatically captured as 
gas internal energy (see Flaig et al.|2010 l. 
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Figure 3. Density-weighted temperature as a function of time. 
Note that the j/-axis has been broken in order to include the 
results from the model at i?, = 1 AU, which has a much higher 
temperature than the other models. 



R = 3 and R — 4 AU, which show a weak trend to cool 
during the whole course of the simulation. 

The time-averaged values of the temperature found in 
in the different simulations are listed in Table [2] Especially 
noticeable is the sharp drop in temperature between the 
simulations at i? = 1 and R = 2 AU, which is due to the 
formation of a dead zone (see below). 



3.1.2 Turbulent activity 

We measure the turbulent activity by calculating the r-cp (or 
x-y) component of the stress tensor normalised to the gas 
pressure, the so-called "alpha parameter", according to the 
following prescription: 

{Txy)p 



(14) 



where T^y = T^y"""^ + T,^"""' and T^y^""^ and T,^"""' are the 
Reynolds and Maxwell stresses, which are defined as 



Rcyn 



■ pVx 5vy and 



-BxBy/flQ (15) 
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Figure 4. Alpha parameter (turbulent stress normalised to gas 
pressure) as a function of time. For further comments, see the 
text. 
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Figure 5. Vertical profiles of the mean magnetic Reynolds num- 
ber, averaged from 60 to 100 orbits. The distance from the mid- 
plane is given in units of the initial pressure scale-height Hq. For 
further comments, see the text. 



where 6vy is the azimuthal component of the gas velocity 
with the velocity of the background shear flow subtracted. 
The quantity Txy is related to the amount of outward angu- 
lar momentum transport taking place at a certain location 
(see, e.g. |Balbus|2003j ). 

The time evolution of the alpha parameter is plotted in 
Fig. |4j the time-averaged values of alpha can be found in 
Table[2] Within the first ten orbits, all simulations reach a 
state of saturated turbulence and remain turbulent until the 
end of the simulation. 

The value of alpha drops sharply when going from R — 
1 AU to 7? = 2 AU, where the ionisation level is very low, 
since the gas is too cool for collisional ionisation and too 
thick for either the the X-rays or the cosmic rays too reach 



Numerical simulations indicate a Schmidt number of order 



the midplane (see Fig. 10 later in the paper). When going 
further outwards, the ionisation level increases, since there 
a significant fraction of the cosmic rays is able to penetrate 
the disc, leading to an increase in the turbulent activity. The 
simulation at i? = 6 AU is again fully turbulent, with a mean 
alpha that is the same as for the simulation at 7? = 1 AU. 
The stresses found in the poorly ionised models at ii = 2 — 4 
AU are one to three orders of magnitude smaller than found 
for the well ionised case. 

Turbulence mixes gas and dust and counteracts the 
force of gravity which causes the dust particles to settle 
towards the midplane. Based on the levels of turbulence 
found in our simulations, we can estimate up to which height 
gas and dust can be considered to be well mixed. We do 
this by equating the settling time-scale = l/Q^t^ (where 
td = /9grCigr/pcs is the drag time) with the turbulent mix- 
ing time-scale tm ~ /Ddust- The dust diffusion coefficient 
f dust is related to the alpha parameter via the Schmidt 



number Sc according to 5c = aCsH / D^ust (Fromang & Pa- 
paloizou]|2006 1. For a density profile in hydrostatic equilib- 
rium, p{z) = poexp{—z^/2H^), the height zq up to which 
gas and dust are well mixed then follows from ta — tm to be 



unity ( Johansen fc Klahr|2005 Turner et al.|20d6 Fromang 
|fc Papaloizou||2006[ ). Setting Sc = 1, and using the initial 
values for the midplane density and sound speed, as well as 
the alpha values from Table [2] Eq. (161 predicts zq — 2.3H 
for the run at R = 2 while for all other runs zq > 3H. 
Since the MRI operates mainly within the first three scale- 
heights (cf. Fig. |6] later in the paper), this means that gas 
and dust can be considered to be well mixed in the MRI 
active regions, so the approximation of a constant dust-to- 
gas ratio is indeed justified. 



3.2 Spatial structure 

3.2.1 Ionisation level & stress profiles 

Concerning the spatial structure of the disc, we first consider 
the question of how the ionisation level and the turbulent 
acitvity vary as a function of position. In Fig. [5] we plot the 
vertical profile of the magnetic Reynolds number, as defined 
in Eq. Q. 

Judging from the condition that regions with Rem > 
can be considered sufficiently ionised for the MRI not to be 
reduced, we can say that the simulations at i? = 1 and 
R — 6 AU can be considered as ideal. All other simulations 
contain at least a small zone of insufficient ionisation near 
the midplane. 

Consequently, the vertical stress profiles, which are 
shown in Fig. [6] are roughly similar for the "ideal" simu- 
lations at i? = 1 and R — 6 AU. The simulations at -R = 5 
and R — 7 AU exhibit only a small dip in the stress near the 
midplane, while in the other simulations, the stresses drop 
noticeably there, forming an extended non-turbulent dead 
zone ranging from a distance of ii = 2 — 4 AU from the 
central star. 



As in the paper of Simon et al. (20111, who also used 



£0 

H 



'2ln—nt4z = 0). 



(16) 



the HLLD Riemann solver, we do not observe a double peak 
in the stress profiles of the well ionised simulations at i? = 
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Figure 6. Vertical profiles of tiie total stress normalised to the 
mean gas pressure {p)p, showing the formation of non-turbulent 
zones near the midplane due to insufficient ionisation. Averages 
have been performed from 60 to 100 orbits. 



1 and R = 6 AU. This is contrary to previous previous 



stratified simulations like that of,Hirose et al. ( 2009 1 or Flaig 



et al. ( 2010 1. Since the only major difference in the numerical 



code as compared to the calculations presented in |Flaig et al.| 
(20101 is the change of Riemann solver, it seems that the 



absence of the double peak profile is indeed due to the use 
of the HLLD solver. 



3.2.2 Structure of the magnetic field 

Concerning the structure of the magnetic field, we plot snap- 
shots of the magnetic field structure for three different runs 
in Fig. [T] The left plot shows the magnetic field structure 
in the fully ideal run at 7? = 1 AU. Within a distance 
of approximately two to three pressure scale-heights from 
the midplane, the MRI is fully operational and the mag- 
netic field is highly tangled. Further outwards, the disc be- 
comes magnetically dominated, which leads to a quenching 
of the MRI, with the magnetic field being predominantly az- 
imuthal. Even further away from the midplane, above four 
pressure-scale heights, the magnetic field becomes distorted 
again. The magnetic field structure that we find for the ideal 
model is very similar to the field structure found byjHirose] 
et al.| ( |2006t (see also [Simon et al.|201lt . 

In the non-ideal runs, the magnetic field structure looks 
quite different, as one can see from the middle plot of Fig.[7j 
which shows the magnetic field structure in the very poorly 
ionised run at J? = 2 AU. Here the magnetic field is laminar 
in the region around the midplane, with only the upper lay- 
ers retaining some level of turbulent activity. This plot can 



(20111 



be compared with Hirose & Turner 

Finally, in the right plot of FigjT we show a snapshot of 
the magnetic field for the run at i? = 5 AU. The magnetic 
field structure is quite similar to that of the run at i? = 1 
AU, with the difference that the magnetic field near the 
midplane is more laminar. 
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Figure 8. Mean vertical temperature profiles, averaged from 60 
to 100 orbits. 
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Figure 9. Mean vertical profiles of the optical depth, averaged 
from 60 to 100 orbits. 



3.2.3 Thermal structure 

Next, we turn to the question of how the thermal structure 
of the disc looks like. We show the temperature profiles in 
Fig. [8] The temperature decreases monotonically with in- 
creasing distance from the central star except at the outer 
edge of the dead zone, where it increases when going from 
J? = 4 AU to i? = 5 AU, due to the increased turbulent 
heating. 

Vertical profiles of the optical depth are shown in Fig.|9] 
In all the simulations, the midplane is optically thick, which 
justifies the use of the one-temperature approximation. Also, 
the photosphere is located well inside the computational do- 
main for all the simulations. 

It is noticeable that the simulations at J? = 2 and R = 3 
AU have a larger central optical depth than the simula- 
tion at i? = 1 AU. Also, the optical depth does not change 
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very much for the five simulations in the range from 2 to 6 
AU. Both observations can be understood from the fact that 
these simulations lie in a temperature range where the opac- 
ity increases strongly with decreasing temperature (scaling 
like K cx T-'^ in this region), which balances the effect of the 
change in surface density on the optical depth. 

3.2.4 Global disc structure 

By interpolating between the results obtained from the local 
simulations at different radii, we can calculate time-averaged 
maps in the R-z plane for the quantities of interest, lead- 
ing to an axisymmetric model of a protoplanetary disc in 
steady state, which is derived from physically realistic three- 
dimensional numerical simulations. 

Fig. |10| shows the total ionisation rate due to the stellar 
X-rays, cosmic rays and the decay of radionuclides. Within 
the first 3-4 AU, the radionuclides dominate the ionisation 
rate near the midplane, while outside of 4 AU, the cosmic 
rays dominate. For the whole radial range considered, the X- 
rays become dominant only in the upper layers, at distances 
> 0.2-0.3 AU away from the midplane. 

Fig.[TT]is intended to give an impression of how the disc 
looks on average in the quasi steady state obtained in the 
local simulations presented in this paper. Shown is the den- 
sity, the location of the photosphere and the magnetosphere 
as well as the location where the fiow becomes supersonic. 
The average extent of the dead zone is also shown, based 
on the criterion that the the magnetic Reynolds number be 
smaller than 1000. The appearance of a second dead zone 
at R = 7 AU is likely an artefact of our particular setup, 
where we neglect the heating due to the stellar irradiation, 
leading to an unrealistically low temperature of only ~ 50 
K and poor ionisation. 



which was located at i? = 1 AU with a surface density about 
six times that of the minimum mass solar nebula, the photo- 
sphere lies always inside the magnetically dominated region, 
although the model presented here is much less massive, and 
we perform simulations at larger radii. As in the jFlaig et aT] 



(20101 model, inside the first 2-3 AU, the turbulence is su- 



personic at the location of the photosphere, but becomes 
slightly subsonic there at larger distances from the central 



star. The conclusion drawn in Flaig et al. (20101 that the 



MRl can be considered as a possible source for the turbulent 
line broadening observed in protoplanetary discs, can thus 
in principle be carried over to models containing dead zones. 
We note that a recent study ( [Simon et al.|[20TT| , also finds 
supersonic velocities above three scale-heights, both for ideal 
models and models containing a dead zone. 

Fig. [12] shows the interpolated mean magnetic field 
strengths found in the simulations presented in this paper. 
The field strengths have of course a minimum at the location 
of the dead zone. Outside of 2 AU, the mean field strengths 
are almost everywhere below ~ 0.2 G. The field strengths 
obtained in the present model are probably too low to ex- 
plain the remnant magnetic fields found in chondritic mete- 
orites, which suggests that a more massive (probably fully 



active) disc is needed to explain these observations ( King & 
Pringle|2010l ). 



We note that as in the model of Flaig et al. (20101, 



4 SUMMARY & CONCLUSION 

We have calculated the spatial structure of a magnetoro- 
tationally turbulent protoplanetary disc from a series of 
physically realistic, local three-dimensional numerical sim- 
ulations. Inside a distance of 1 AU from the central star, 
the disc gas is sufficiently ionised by collisional ionisation, 
so the disc midplane is fully turbulent. Between 1 and 2 AU, 
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Figure 10. lonisation rate due to the external ionisation sources. Tlie black line corresponds to the boundary between regions where 
diflferent ionisation sources dominate. 
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Figure 11. Global spatial structure. Shown is a contour plot of the density. The shaded area corresponds to the region where the 
magnetic Reynolds number is smaller than 1000. The red line marks the location of the photosphere, the black line corresponds to the 
the location where magnetic pressure and gas pressure are equal and the green line corresponds to the boundary where the flow starts 
to become supersonic. 



the temperature drops below the threshold of 900 K which 
is needed for collisional ionisation to be effective, leading 
to a transition from the fully turbulent state to a state that 
contains a central non-turbulent dead-zone. When going fur- 
ther outwards, the disc gas becomes thinner, which allows 
a larger fraction of the stellar X-rays to reach the region 
around the midplane, leading to a gradual increase in the 
ionisation level and also the turbulent activity. At i? = 6 
AU, the midplane is again fully turbulent. 

The disc gas is heated internally (via turbulent and 
Ohmic dissipation) with the heating due to the stellar ir- 
radiation being neglected. The simulations reach (approxi- 



mate) thermodynamical equilibrium, except possibly for the 
simulation located at 7? = 3 and R — 4 AU, which continue 
to cool slowly for the whole course of the simulation. For 
these two particular runs, the temperature may therefore 
dependent on the initial conditions, while for all the other 
simulations, the temperature comes out self-consistently as 
a result of the balance between internal heating and radia- 
tive cooling. Since in most of our simulations the disc gas is 
for a large part optically thick, we do not expect the stellar 
irradiation to change the temperature in the interior of the 
disc drastically, except for the simulation at R = 7, where 
the temperature is indeed unrealistically low (when com- 
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Figure 12. Mean magnetic field strength. 



pared, for example, with the passive disc model of |Chiangl 
fc Goldreichfl997| . 

It should be noted that the use of the local approxi- 
mation is a major restriction of the present work, since lo- 
cal simulations cannot capture the full dynamics of a global 
simulation (see, for example, Sorathia et al.|2011 1. However, 
apart from the fact that in a global setup the surface den- 
sity would change over time, the spatial structure in a global 
model using the same physics will likely not be dramatically 
different from the structure coming out in the local frame- 
work. Therefore, the approach taken in the present work is 
a useful first step in deriving the structure of protoplanetary 
discs from physically realistic numerical simulations. 

To summarise, the simulations presented in this work 
are able to capture all the basic physics that is thought to 
be important for protoplanetary disc dynamics, and do not 
suffer from the lack of self-consistency of previous isothermal 
models. For the first time, it is therefore possible to obtain 
a picture of the global structure of a magnetorotationally 
turbulent protoplanetary disc from first-principles numerical 
simulations. 
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